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Abstract 

We are concerned with the numerical solution obtained by splitting methods of certain parabolic partial 
differential equations. Splitting schemes of order higher than two with real coefficients necessarily involve 
negative coefficients. It has been demonstrated that this second-order barrier can be overcome by using 
splitting methods with complex-valued coefficients (with positive real parts). In this way, methods of 
orders 3 to 14 by using the Suzuki- Yoshida triple (and quadruple) jump composition procedure have 
been explicitly built. Here we reconsider this technique and show that it is inherently bounded to order 14 
and clearly sub-optimal with respect to error constants. As an alternative, we solve directly the algebraic 
equations arising from the order conditions and construct methods of orders 6 and 8 that are the most 
accurate ones available at present time, even when low accuracies are desired. We also show that, in the 
general case, 14 is not an order barrier for splitting methods with complex coefficients with positive real 
part by building explicitly a method of order 16 as a composition of methods of order 8. 

Keywords: composition methods, splitting methods, complex coefficients, parabolic evolution equations. 
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1 Introduction 

In this paper, we consider linear evolution equations of the form 

dix 

— {t) = Au{t)+Bu{t), u{Q) = uo, (1.1) 
dt 

where the (possibly unbounded) operators A , B and A + B generate semi-groups for positive t over 
a finite or infinite Banach space X . Equations of this form are encountered in the context of parabolic 
partial differential equations, a prominent example being the inhomogeneous heat equation 

du 

-^(x, t) = Au{x, t) + V{x)u{x, t), 
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where t > , x G M'^ or x G T'^ and A denotes the Laplacian in x . 

A method of choice for solving numerically (1.1) consists in advancing the solution alternatively along 
the exact (or numerical) solutions of the two problems 

^(t) = Au{t) and ^{t) = Bu{t). 

Upon using an appropriate sequence of steps, high-order approximations can be obtained -for instance with 
exact flows- as 

i^[h}) = Q^^oB gfeaiA g/ifeiB . . . ^hasA ^hb^B q 2) 

The simplest example within this class is the Lie-Trotter splitting 

^hA^hB ehB^hA^ (J 3) 

which is a first order approximation to the solution of ( 1 . 1 ), while the symmetrized version 

5(/i) =e'^/2^e'^^e'^/2^ or = e'^/^B ^hA ^^2^ (14) 

is referred to as Strang splitting and is an approximation of order 2 . 

The application of splitting methods to evolutionary partial differential equations of parabolic or mixed 
hyperbolic-parabolic type constitutes a very active field of research. For this class of problems it makes 
sense to split the spatial differential operator, each part corresponding to a different physical contribution 
(e.g., reaction and diffusion). Although the formal analysis of splitting methods in this setting can be car- 
ried out by power series expansions (as in the case of ordinary differential equations), several fundamental 
difficulties arise, however, when establishing convergence and rigorous error bounds for unbounded opera- 
tors [HKLRIO]. Partial results exist for hyperbolic problems [TT95, Tang98, HKLRIO], parabolic problems 
[DS02, HV03] and for the Schrodinger equation [JLOO, Lub08], even for high order splitting methods [Tal08]. 

In [HO09a], it has been established that, under the two conditions stated below, a splitting method of the 
form (1.2) is of order p for problem (1.1) if and only if it is of order p for ordinary differential equations in 
finite dimension. In other words, if and only if the difference ^'(/i) — e'^'^^"'"^) admits a formal expansion of 
the form 

^{h) - e^i^+B) = hP+^Ep+i + + ■■■ (1-5) 

The two referred conditions write (see [HO09a] for a complete exposition): 

1. Semi-group property: A, B and A-\-B generate semi-groups on X and, for all positive t, 

|[e*^|[ < e'^^*, ||e*^|| < e'"^* and ||e*(^+^)|| < e"* 
for some positive constants ua , w b and u) . 

2. Smoothness property: For any pair of multi-indices (ii, . . . , im) and (ji, . . . ,jm) with ii + • • • + 

«m + ii H h jm = P + 1 , and for all t G [0, T] , 

ll^'i^-^'i . . . A'^B^"' e*(^+^)no|| < C 

for a positive constant C . 

However, designing high-order splitting methods for (1.1) is not as straightforward as it might seem at first 
glance. As a matter of fact, the operators A and B are only assumed to generate C° semi-groups (and not 
groups). This means in particular that the flows e*"^ and/or e*^ may not be defined for negative times (this 
is indeed the case, for instance, for the Laplacian operator) and this prevents the use of methods which embed 
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negative coefficients. Given that splitting methods with real coefficients must have some of their coefficients 
ai and bi negative' to achieve order 3 or more, this seems to indicate, as it has been believed for a long 
time within the numerical analysis community, that it is only possible to apply exponential splitting methods 
of at most order p = 2 . In order to circumvent this order-barrier, the papers [HO09b] and [CCDV09] si- 
multaneously introduced complex-valued coefficients^ with positive real parts. It can indeed be checked in 
many situations that the propagators e^^ and e^^ are still well-defined in a reasonable distribution sense 
for z G C, provided that > . Using this extension from the real hne to the complex plane, the 

authors of [HO09b] and [CCDV09] built up methods of orders 3 to 14 by considering a technique known 
as triple-jump composition^ and made popular by a series of authors: Creutz & Gocksch [CG89], Forest 
[For89], Suzuki [Suz90] and Yoshida [Yos90]. 

In this work, we continue the search for new methods with complex coefficients with positive real parts. 
Eventually, our objective is to show that, compared to the methods built in [HO09b] and [CCDV09] by 
applying the triple-jump (or quadruple-jump) procedure, it is possible to construct more efficient methods 
and also of higher order by solving directly the polynomial equations arising from the order conditions. In 
particular, we construct methods of order 6 and 8 with minimal local error constants among the methods 
with minimal number of stages. We also construct a method of order 16 , obtained as a composition based 
on an appropriate 8th order splitting method. 

Here we are particularly interested in obtaining new splitting methods for reaction-diffusion equations. 
These constitute mathematical models that describe how the population of one or several species distributed 
in space evolves under the action of two concurrent phenomena: reaction between species in which predators 
eat preys and diffusion, which makes the species to spread out in space"^. From a mathematical point of view, 
they belong to the class of semi-linear parabolic partial differential equations and can be represented in the 
general form 

— = DAn + F(u), 

where each component of the vector u{x, t) € M'^ represents the population of one species, D is the real 
diagonal matrix of diffusion coefficients and F accounts for all local interactions between species. Strictly 
speaking, the theoretical framework introduced in [HO09b] does not cover this situation if F is nonlinear, 
so that (apart from section 4, where we successfully integrate numerically an example with nonlinear F ) we 
will think of F as being linear. The important feature of ^ = DA here is that it has a real spectrum: hence, 
any splitting method involving complex steps with positive real part is suitable for that class of problems. In 
principle, one could even consider splitting methods with 's having positive real part and unconstrained 
complex bi 's. 

It may be worth mentioning that the size of the arguments of the Oj coefficients of the splitting method 
is a critical factor when the diffusion operator involves a complex number, for instance, an equation of the 
form 

f)ii 

— = §Au + Fiu), (1.6) 

where 5 is a complex number with ^{5) > . The choice F{u) = fi^u^ + /X2U^ + jiiu + is known 
as the cubic Ginzburg-Landau equation [FT88]. In this situation, the values of the := arg((5) + arg(aj) 
determine whether the splitting method makes sense for this specific value of 5 . If for all i = 1, . . . , s , 

'The existence of at least one negative coefficient was siiown in [Slie89, Suz91], and the existence of a negative coefficient for 
both operators was proved in [GK96]. An elementary proof can be found in [BC05]. 

^Methods with complex-values coefficients have also been used in a similar context [Ros63] or in celestial mechanics [Cha03]. 
And its generalization to quadruple-jump. 

''Apart from biology and ecology, systems of this sort also appear in chemistry (hence the term reaction), geology and physics. 

^The choice F{u) — u{l — u) yields Fisher's equation and is used to describe the spreading of biological populations; the 
choice F[u) — u{l — v?) describes Rayleigh-Benard convection; the choice F{u) — u{l — u)[u — a) with < a < 1 arises 
in combustion theory and is referred to as Zeldovich equation. 
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fflj € [— f , +f ] then the method is well defined. In order for the method to be applicable to such class of 
equations, it would make sense trying to minimize the value of maxj=i ^ | arg(aj)| . 

In this work, however, we focus on the case where the operator A has a real spectrum, and thus we will 
only require that | arg(aj)| < 7r/2 (i.e., K(aj) > ) while minimizing the local error coefficients. In this 
sense, we have observed that the coefficients of accurate splitting methods with 5?(aj) > tend to have also 
bi coefficients with positive real parts. 

The plan of the paper is the following. In Section 2, we shall prove that if an s -jump construction is 
carried out from a basic symmetric second-order method, it is bounded to order 14 and no more and we 
will further justify why solving directly the system of order conditions leads to more efficient methods. In 
Section 3, we solve the corresponding order conditions of methods based on compositions of second order 
integrators and construct several splitting methods whose coefficients have positive real part. In particular, 
in subsection 3.3 we present splitting methods of orders 6 and 8 , obtained as a composition scheme with 
Strang splitting as basic integrator. In addition, with the aim of showing that 14 is not an order barrier in 
general, we have built explicitly a method of order 16 as a composition of methods of order 8, which in 
turn is obtained by composing the second order Strang splitting. In section 4 we describe the implementation 
of the various methods obtained in this paper and show their efficiency as compared to already available 
methods on two test problems. Finally, section 5 contains some discussion and concluding remarks. 

2 An order barrier for the s -jump construction 

A simple and very fruitful technique to build high-order methods is to consider compositions of low-order 
ones with fractional time steps. In this way, numerical integrators of arbitrarily high order can be obtained. 
For sphtting methods aimed to integrate problems of the form (1.1), it is necessary, however, that the coeffi- 
cients have positive real part. The procedure has been carried out in [CCDV09, HO09b], where composition 
methods up to order 14 have been constructed. We shall prove here that 14 constitutes indeed an order barrier 
for this kind of approach. In other words, the composition technique used in [CCDV09, HO09b] does not 
allow for the construction of methods having all their coefficients in C+ := {z € C : $R(z) > 0} with orders 
strictly greater than 14 . With this goal in mind we consider the following two families of methods: 



Family I. Given a method of order p, $ W (h) = + 0{hP+^) , a sequence of methods of orders 

p + l,p + 2, . . . can be defined recursively by the compositions 

<I>[f+5l(/,)= JJ$b+'?-i](a^^^/,) .= $b+a-i](c,^^^;,) ... ^\p+i-n(^a,,^h), (2.1) 
i=l 

where for all q>l, 

rriq ruq 

(VI <i<mg, a,,, 7^0), "^'^ = ^ Z]<? = 0- 

i=l 1=1 

(Hereafter, we will interpret the product symbol from left to right). Notice that if p + g is even, the second 
condition has only complex solutions. 



Family II. Given a symmetric method of order 2p , <|[^p] (h) , a sequence of methods of orders 2{p+l) , 
2{p + 2), . . . can be defined recursively by the symmetric compositions 

niq 

V g > 1, |.[2(P+'?)1 (/i) = Jl $[2(p+g)-2] ./i) (2.2) 

i=l 
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where a<j,mg+i-i = ag,i, i = 1, 2, . . . , and for all g > 1 , 

rUq 1 

( V 1 < z < mq, a,,, / 0) , Yl = 1 ^"^d "'ir*^^' = 0. 



,2{p+g)+i 

c 

=1 i=l 



Methods of this class with real coefficients have been constructed by Creutz & Gocksch [CG89], Suzuki 
[Suz90] and Yoshida [Yos90]. However, the second condition clearly indicates that at least one of the coeffi- 
cients must be negative. In contrast, there exist many complex solutions with coefficients in . 
Generally speaking, starting from #[^1 (h) , the {p + 1) -th member of family I is of the form 

rrip / mp-i / / "ii \ \ \ 

$^+i](/i)=j] n ... n^^''K..«,-i,,_, ... (2.3) 

ip=l \ip_i=l V \ii=l / / / 

and has coefficients 

p 

JJoj-.i, I < ii <mi, . . . ,1 < ip < lUp. (2.4) 
i=i 

A similar expression holds, of course, for methods of family II, starting from l>[^l(/i) . Symmetric 
compositions for the cases mi = m2 = ■ ■ ■ = nip = 3 and mi = m2 = ■ ■ ■ = nip = 4 , correspond to the 
triple and quadruple jump techniques, respectively. 

Lemma 2.1 Let ^{h) be a consistent method (i.e., a method of order p > 1 ) and assume that the method 

I 

^(/i) = J]$(Qi/i) (2.5) 

1=1 

is also consistent (i.e., = 1)- If there exists k, 1 < k < I , such that 5ft(afc) < 0, then any consistent 

method of the form 

m m / I \ 

ll^{f3jh) = llill^f^,aih)\ (2.6) 

j=i j=i \i=i J 

has at least one coefficient fijUk , ^ 1^ j < tti , such that ^{(3jak) < . 
Proof: By consistency one has X^JLi /3j = 1 , so that 

m / in \ / m \ 

j=i \i=i / V i=i / 



This implies the statement. □ 

Lemma 2.2 For k > 2 and 

i = 1, . . . ,k . Then we have 



Lemma 2.2 For k > 2 and r > 2 , consider (zi, . . . , z^) € (C+)'^ such that "^^^i zl = 0, Zi ^ for 

TT 

max Argizi) — min Arg(2;j) > — . 
i=l,...,fc i=i,...,fc r 



Proof: All the Zi 's belong to the sector K„{9) = {z e C : \a - Arg(z)| < 9} with 

a = — ( max Arg(zj) + min Arg(2;j) ] and ^ — ~ ( max Arg(zj) — min Arg(2;j 

where Arg is the principal value of the argument (see the left picture of Figure 1). Now, assume that 9 < ^ . 
Then, all the 's belong to Kra{r9) , which, since < ^ , is a convex set. This implies that Yli=i 
also belongs to Kra{r9) (see the right picture of Figure 1). The inequality r9 < ^ being strict and the Zi 's 
non-zero, we have furthermore Yli=i ¥^ ^ ^ which contradicts the assumption. The result follows. □ 
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Figure 1: The enveloping sectors of {zi, . . . , Zk} C C-|- and of {z^, . . . , z^.} C C (for r = 2 ). 

Theorem 2.3 (i) Starting from a first-order method <^^^\h) , all methods <^^^^^\h) of order p+1, p = 
3, 4, . . . from Family I have at least one coefficient with negative real part, (ii) Starting from a second-order 
symmetric method $^^1 (h) , all methods |>[^p+^] (/i) of order 2p + 2, p = 7,8, . . . from Family II have at 
least one coefficient with negative real part. 

Proof: We prove at once the two assertions. We first notice that, according to Lemma 2.1, if method <^^P\h) 
of Family I (respectively, method l>[^^](/i) of Family II), has a coefficient with negative real part, then all 
subsequent methods , (7 > 1 , of Family I (respectively, methods <|[2(p+'?)] (/i) of Family II), also 

have a coefficient with negative real part. Hence, we assume that all methods <^^'^~^^\h) , q = 1, . . . ,p 
from Family I (respectively, all methods of Family II), have all their coefficients in C+ . Using 

Lemma 2.2 we have 



Mq=l,...,p, max Arg(Q;g,i) 

1=1,. ..,mq 



TT 



mill Arg(ag,i) > 

1=1,. ..,mq Q ~r i- 



(respectively 



so that 



yq = l,...,p, max Arg(a, 



TT 



min Arg(aq,i) > 

1=1, ...,mq 2q -\- L 



max Arg IT aj^., 

\J=1 



mm 

l<ii<mi,...,l<jp<m, 



Arg I Yl 



TT 

I > - + 



TT 



P+1 



(respectively 



max Arg ( TT a,- j . | — min Arg ( TT a, i . 1 > — h • • • H )• 

l<il<mi,...,l<ip<mp \ -l I ^<ii<m,i,...,l<ip<mp \ ~q / 3 2^+1 



Now, since | + | + |>l,p = 3in the first case and thus the first statement follows. For Family II, since 



i + i + -- - + -rr>l, then p = 7 , thus leading to the second statement. □ 



Remark 2.4 No method of Family I with coefficients in C-(_ can have an order strictly greater than 3. Such 
methods of order 3 have been constructed in [HO09b]. Similarly, no method of Family II with coefficients in 
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C+ can have an order strictly greater than 14 . Such methods with orders up to 14 have been constructed 
in [CCDV09, HO09b ]. For instance, let us consider the quadruple jump composition 

$[2p+2](^) = ^i^P^ap^^h) <|[2p](ap 2/j) |.[2p](ap 2/i) 4>[2p] (ap,i/i), (2.7) 

where ap^i + ap_2 = 1/2 and + o?p^2'^ = 0. This system has p solutions (and their complex 

conjugate). The solution with minimal argument is, as noticed in [CCDV09, HO09b], 




ap,i = - i + t p r- , ap,2 = ap,i 



(and its complex conjugate). It is straightforward to verify that arg(ap i) = 2(2p+i) ''"'^ arg(ap^2) 

-2(2^'^«^^«^ 

arg(ap,i) - arg(ap,2) 



2p+l 

Comparing with the proof of Theorem 2.3, we observe that the bounds obtained there are sharp since a family 
of methods do exist satisfying the equality. 

Remark 2.5 It is however possible to construct a composition method with all coefficients having positive 
real part of order strictly greater than 14 directly from a symmetric second order method. For example, in 
Subsection 3.3 we present a new method of sixteenth-order built as 

21 15 

<^>^^^\h) = Yl^^^^{aih), with ^^^^{h) = ll^^^^{^jh) (2.8) 

1=1 j=i 

and the coefficients satisfying ^{ai(3j) > for all i = 1,...,21, j = 1,...,15, with 022-1 = 
ttj, /3i6-j = /3j, hj = 1,2,.... Here, ^^^\h) is a symmetric composition of symmetric second order 
methods, but it is not a composition of methods of order 4 or 6, and similarly for <b^'^^\h) , which is not a 
composition of methods of orders 10, 12 or 14. 

Theorem 2.6 Splitting methods of the class (1.2) with K(ai) > exist at least up to order 44. 

Proof: In [CCDV09], a fourth-order method was obtained with Oj € M+ . In a similar way, we have also 
built a sixth-order scheme with Oi G (whose coefficients can be found at 

http : / /www .gicas. uji.es/ Re search / splitting-complex . html). Using this as the basic 
method in the quadruple jump (2.7) with coefficients chosen with the minimal argument and, since 

1111 

^ + • • • + 7T7 < 1< 7; + • • • + 77 

7 43 7 45 

we conclude that all methods obtained up to order 44 will satisfy that K(aj) > . □ 

Remark 2.7 The question of the existence of splitting schemes at any order with K(ai) > remains still 
open. 



3 Splitting methods with all coefficients having positive real part 
3.1 Order conditions and leading terms of local error 

We have seen that the composition technique to construct high order methods inevitably leads to an order 
barrier. In addition, the resulting methods require a large number of evaluations (i.e. 3""^ or 4"^^ eval- 
uations to get order 2n using the triple or quadruple jump, respectively) and usually have large truncation 
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errors. In this section we show that, as with real coefficients, it is indeed possible to build very efficient high 
order splitting methods whose coefficients have positive real part by solving directly the order conditions 
necessary to achieve a given order p . These are, roughly speaking, large systems of polynomial equations 
in the coefficients aj , bi of the method (1.2), arising by requiring that the formal expansion of the method 
satisfies (1.5) for arbitrary non-commuting operators A and B. 

Different (but equivalent) formulations of such order conditions exist in the literature [HLW06]. Among 
them, the one using Lyndon multi-indices is particularly appealing. It was first introduced in [CM09] (see 
also [BCM08]) and can be considered as a variant of the classical treatment obtained in [MSS99]. 

This analysis shows that the number of order conditions for general splitting methods of the form (1.2) 
grows very rapidly with the order p , even when one considers only symmetric methods. For instance, 
there are 26 independent 8th-order conditions and 82 lOth-order conditions for a consistent symmetric 
splitting method. It makes sense, then, to examine alternatives to achieve order higher than six. This can 
be accomplished by taking compositions of a basic symmetric method of even order In particular, if we 
consider any of the two versions of Strang splitting (1.4) as the basic method S{h) , then, for each 7 = 
(7i,...,7m) eC-, 

^{h) = S{^ih)---S{jsh) (3.1) 

will be a new splitting method of the form (1.2). Now the consistency condition reads 

71 + • • • + 7s = 1- (3.2) 

As for the additional conditions to attain order p , these can be obtained by generalizing the treatment done in 
[MSS99]. Splitting methods with very high order can be constructed as (3.1) by considering as basic method 
S{h) a symmetric method of even order 2g > 2 . In that case, it can be shown that the corresponding number 
of order conditions is considerably reduced with respect to (1.2). Thus, for instance, if S{h) is a symmetric 
splitting method of order eight, a method (3.1) satisfying the consistency condition (3.2) and the symmetry 
condition 

7s-j+i=7j, l<j<s (3.3) 
only needs to satisfy 10 additional conditions to attain order sixteen. 

3.2 Leading term of the local error 

To construct splitting methods of a given order p within a family of schemes, we choose the number s 
of stages in such a way that the number of unknowns equals the number of order conditions, so that one 
typically has a finite number of isolated (real or complex) solutions, each of them leading to a different 
splitting method. Among them, we will be interested in methods such that, either ai > (and each bi 
are arbitrary complex numbers) or 5?(aj) > and > . The relevant question at this point is how 

to choose the 'best' solution in the set of all solutions satisfying the required conditions. It is generally 
accepted that good splitting methods must have small coefficients ai,bi . Methods with large coefficients 
tend to show bad performance in general, which is particularly true when relatively long time-steps h are 
used. In addition, as mentioned in the introduction, when applying splitting methods to the class of problems 
considered here, the arguments of the complex coefficients Oj, bi must also be taken into account. 

In order to chose the best scheme among two methods with coefficients of similar size included in sectors 
with similar angle, we analyze the leading term of the local error of the splitting method. If (1.2) is of order 
p , then we formally have that 

Ep+i := J2 v^,_^,^{-i)A'^B^^...A'^--^B^^-, 

ilH |-j2m=P+l 
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where 7 = (71, . . . , 7^) is given in terms of the original coefficients Oj , bi of the integrator by 



7i-i + 7j 

«i=7j, ftj-i = ^^ 

(70 = 7s+i = ) and each Wii,...,i2„(7) is a linear- combination of polynomials in 7 [BCM08]. 

Incorporating that into the results in [HO09b], it can be shown that, if the smoothness assumption stated 
in the introduction is increased from p+ 1 to p+2 , then for sufficiently small h , the local error is dominated 

by ||/if+ii?p+ino||- 

Our strategy to select a suitable method among all possible choices is then the following: first choose a 
subset of solutions with reasonably small maximum norm of the coefficient vector (71, . . . , 7^) , and then, 
among them, choose the one that minimizes the norm 

E 1^^1-^2^(7)1 (3-4) 

iiH \-i2m=P+l 

of the coefficients of the leading term of the local error. This seems reasonable if one is interested in choosing 
a splitting method that works fine for arbitrary operators A and B satisfying the semi-group and smoothness 
conditions mentioned in the introduction. Of course, this does not guarantee that a method with a smaller 
value of (3.4) will be more precise for any A and B than another method with a larger value of (3.4). 
However, we have observed in practice when solving the order conditions of different families of splitting 
methods, that the solution that minimizes (3.4) tend to have smaller values of most (or even all) coefficients 
bn - i2m(7)l when compared to a solution having a larger norm (3.4) of the coefficients of the leading term 
of the local error. 

When A and B are operators in a real Banach space X , then it makes sense to compute the approxi- 
mations u„ = u{tn) , as Un = . In that case, the argument above holds with ^(/i) replaced 
by ^{h) = and the local error coefficients Wji,.. .,12^(7) replaced by 5?(t'ii,...,j2m(7)) • In that 
case, (3.4) should be replaced by 

E I^K,...,i2™(7))l (3.5) 

iiH \-i2m=P+l 

as a general measure of leading term of the local error. 

3.3 High order splitting methods obtained as a composition of simpler methods 

Order 6. We first consider sixth-order symmetric splitting methods obtained as a composition (3.1) of the 
Strang splitting (1.4) as basic method. In this case the coefficients 7^ must satisfy three order conditions, in 
addition to the symmetry (3.3) and consistency requirements, to achieve order six. We thus take s = 7 , so 
that we have three equations and three unknowns. Such a system of polynomial equations has 39 solutions in 
the complex domain (three real solutions among them), 12 of them giving a splitting method with coefficients 
of positive real part. According to the criteria described in Subsection 3.2, we arrive at the scheme 

71 = 77 = 0.116900037554661284389 +0.043428254616060341762i (3.6) 

72 = 76 = 0.12955910128208826275 - 0.12398961218809259330i, 

73 = 75 = 0.18653249281213381780 + 0.00310743071007267534i, 

74 = 0.134016736702233270122 + 0.1549078537239191523961 

This method turns out to correspond to one of those obtained by Chambers (see Table 4 in [Cha03]). 

Order 8. For consistent symmetric methods (3.1) of order eight, we have seven order conditions. By 
taking s = 15 stages, one ends up with a system of seven polynomial equations and seven unknowns. We 
have performed an extensive numerical search of solutions with small norm, finding 326 complex solutions. 
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Among them, 162 lead to splitting methods whose coefficients possess positive real part. The best method, 
according to the criteria established in Subsection 3.2, is 



71 = 


715 


= 0, 


,053475778387618596606 


+ 0.006169356340079532510i 


72 = 


714 


= 0, 


,041276342845804256647 


- 0.06994857439070781495H 


73 = 


713 


= 0, 


,086533558604675710289 


- 0.023112501636914874384^ 


74 = 


712 


= 0, 


.079648855663021043369 


+ 0.049780495455654338124Z 


75 = 


711 


= 0, 


,069981052846323122899 


- 0.052623937841590541286i 


76 = 


710 


= 0, 


,087295480759955219242 


+ 0.010035268644688733950^ 


77 = 


= 79 


= 0, 


,042812886419632082126 


+ 0.076059456458843523862^ 




78 


= 0, 


,077952088945939937643 


+ 0.007280873939894204350Z 



Order 16. Motivated by the results in Section 2, we have also constructed a splitting method of order 16. 
Our aim, rather than proposing a very efficient scheme, is to show that the barrier of order 14 existing for 
methods built by applying the recursive composition technique starting from order two (family II) does not 
apply in general. 

The construction procedure can be summarized as follows. We consider a consistent symmetric com- 
position of the form (3.1), where now the basic method S{h) is any symmetric eighth order scheme. 
Under such conditions, ten order conditions must be satisfied to achieve order 16. We accordingly take 
s = 21 , so that we have ten polynomial equations with ten unknowns. We have performed an extensive 
numerical search of solutions with relatively small norm, finding 70 complex solutions with positive real 
part. Combined with the 162 methods of order eight, this leads to 11340 different 16 -th order splitting 
methods with s = 315 stages. Among them, only 324 give rise to splitting methods with coefficients 
of positive real part. The coefficients of the method that we have determined as optimal can be found at 
http : / /www .gicas.uji.es/ Re search/ splitting-complex . html. 

Remark 3.1 Notice that, whereas for the sixth-order method (3.6) one has 

vr vr 

max Arff(7j) — min Are;(7j) ~ 1.621 < 1 

i=l,...,7 i=l,...,7 3 5 

and then the coefficients 7^ are distributed in a narrower sector than for triple or quadruple jump methods, 
for the eighth-order method ( 3. 7) one has 

vrTTvr 

max Arg(7j) — min Arg(7j) ~ 2.5997 > 1 1 . 

i=l,,,,,15 i=l,...,15 3 5 7 

This method, whereas being the most efficient, is not the appropriate one to be used as basic scheme to 
build higher order methods by composition. The previous 16th-order integrator has been built starting with 
another 8th-order method whose coefficients are placed in a narrower sector 



4 Numerical tests 

For our numerical experiments, we consider two different test problems: a linear reaction-diffusion equation, 
and the semi-linear Fisher's equation, both with periodic boundary conditions in space. It should be men- 
tioned here that this last case is not covered by the theoretical framework summarized in the Introduction. For 
each case, we detail the experimental setting and collect the results achieved by the different schemes. Our 
main purpose here is just to illustrate the performance of the new splitting methods to carry out the time inte- 
gration as compared with those constructed by using the Yoshida-Suzuki triple jump composition technique 
for both examples. Notice that, in this sense, the particular scheme used to discretize in space is irrelevant. 
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For that reason, and to keep the treatment as simple as possible, we have applied a simple second-order finite 
difference scheme in space. 

The numerical approximations u„ obtained by each method ^'(/i) are computed as n„ = . 
In other words, we project on the real axis after completing each time step. 



4.1 A linear parabolic equation 

Our first test-problem is the scalar equation in one-dimension 

du(x, t ' 



dt 



aAu{x,t) + V{x,t)u{x,t), u{x , 0) = uo{x) , (4.1) 



with uo{x) = sin(27rx) and periodic boundary conditions in the space domain [0, 1] . We take a = j , 
V{x, t) = 3 + sin(27ra;) and discretize in space 

Xj=j{6x), j = l,...,N with 6x = l/N, 

thus arriving at the differential equation 

^ = aAU + BU, (4.2) 
dt 

where U = (ni, . . . ,un) G . The Laplacian A has been approximated by the matrix A of size N x N 
given by 

/-2 1 1 \ 

1 -2 1 

A = ^ 1 -2 1 

{Sx)^ 

\l 1-2/ 

and B = diag(F(xi), . . . , V{xn)) . We take = 100 points and compare different composition methods 
by computing the corresponding approximate solution on the time interval [0,1] . In particular, we consider 
the following schemes: 

• Strang: The second-order symmetric Sti"ang splitting method (1.4); 

• (TJ6): The sixth-order triple jump method (Proposition 2.1 in [CCDV09]) based on Strang's second- 
order method; 

• (TJ6A): The sixth-order triple jump method (Proposition 2.2 in [CCDV09]) based on Strang's second- 
order method; 

• (TJ8A): The eighth-order triple jump method (Proposition 2.2 in [CCDV09]) based on Strang's second- 
order method; 

• (P6S7): The sixth-order method (3.6); 

• (P8S15): The eighth-order method (3.7). 

We compute the error of the numerical solution at time t = 1 (in the 2 -norm) as a function of the number of 
evaluations of the basic method (the Strang splitting) and represent the outcome in Figure 2. In the left panel 
we collect the results achieved by the Strang splitting and the previous sixth-order composition methods, 
whereas the right panel corresponds to eighth-order methods. We have also included, for reference, the curve 
obtained by (P6S7). 

The relative cost (w.r.t. Strang) of a method composed of s steps is approximated by 4s , where the 
factor 4 stands here for an average ratio between the cost of complex arithmetic compared to real arithmetic. 
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A remarkable outcome of these experiment is that methods (P6S7) and (P8S15) outperform Strang's splitting 
(and actually all other methods tested here) even for low tolerances. Scheme (P8S15), in particular, proves 
to be the most efficient in the whole range explored. The gain with respect to triple jump methods is also 
significant and completely support the approach followed here. 




Figure 2: Error versus number of steps for the linear reaction-diffusion equation (4. 1). 



4.2 The semi-linear reaction-diffusion equation of Fisher 

Our second test-problem is the scalar equation in one-dimension 

du{x,t) ^ ^u{x,t)+F{u{x,t)), u{x,0) = uo{x), (4.3) 
with periodic boundary conditions in the space domain [0, 1] . We take, in particular, Fisher's potential 

F{u) =u{l-u). 

The splitting considered here corresponds to solving, on the one hand, the linear equation with the operator 
A being the Laplacian, and on the other hand, the nonlinear ordinary differential equation 

— = u{x, t)(l - u{x, t)) 

with initial condition 

u{x, 0) = uo{x). 

Note that it can be solved analytically as 

(e* - 1) 



u{x, t)) = uo{x) + uo{x){l - uoix)) 



l + no(x)(e* - 1)^ 



which is well defined for small complex time t . We proceed in the same way as for the previous linear case, 
starting with uo(x) = sin(27rx) . After discretization in space, we arrive at the differential equation 

^=AU + F{U), (4.4) 
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where [/ = (ui, . . . , uat) e and F{U) is now defined by 

F{U) = (ni(l - ui), un{1 - un)) ■ 

We clioose N = 100 and compute the error (in the 2 -norm) at the final time t = 1 by applying the same 
composition methods as in the linear case. The results are collected in Figure 3, where identical notation has 
been used. Notice that, strictly speaking, the theoretical framework upon which our strategy is based does 
not cover this nonlinear problem. Nevertheless, the results achieved are largely similar to the linear case. In 
particular, the new 8th-order composition method is the most efficient even for moderate tolerances. 

Error versus cost in log-log scale Error versus cost in log-log scale 




2.4 2.6 2.8 3 3.2 3.4 3.6 2.4 2.6 2.8 3 3.2 3.4 3.6 



Number of steps of ttie Strang splitting mettiod Number of steps of ttie Strang splitting mettiod 



Figure 3: Error versus number of steps for the semi-linear reaction-diffusion equation (4.3). 



5 Concluding remarks 

Splitting methods with real coefficients for the numerical integration of differential equations of order higher 
than two have necessarily some negative coefficients. This feature does not suppose any special impediment 
when the differential equation evolves in a group, but may be unacceptable when it is defined in a semi-group, 
as is the case with the evolution partial differential equations considered in this paper. One way to get around 
this fundamental difficulty is to consider splitting schemes with complex coefficients having positive real part. 
This has been recently proposed for diffusion equations in [CCDV09, HO09b]. Splitting and composition 
methods with complex coefficients have been considered in different contexts in the literature (see [BCMIO] 
and references therein). 

In [CCDV09, HO09b], splitting methods up to order 14 with complex coefficients with non-negative real 
part have been recursively constructed either by the so-called triple-jump compositions or by the quadruple- 
jump compositions, starting from the symmetric second-order Strang splitting. In this work we prove that 
there exists indeed an order barrier of 14 for methods constructed in this way. More generally, we show 
that no method of order higher than 14 with coefficients having non-negative real part can be constructed 
by sequential s -jump compositions starting from a symmetric method of order 2. We further show, by 
explicitly obtaining methods of order 16 (as the composition of a basic symmetric method of order 8), that 
this order barrier does not apply for general composition methods (non-necessarily constructed by recursive 
applications of s -jump compositions) with complex coefficients with non-negative real part. 

In addition to this order barrier, another drawback of methods resulting from applying the s -jump com- 
position procedure is that for high orders they require a larger number of stages (i.e. number of compositions 
of the basic symmetric second order method) than methods obtained by directly solving the order conditions 
with the minimal number of stages. For instance, methods of order 6 (respectively, 8) obtained with triple 
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jump compositions need 9 (resp. 27) compositions of the basic second order method, whereas, as we show 
in the present work, methods of order 6 (resp. 8) can be constructed (by directly solving for the required 
order conditions) with 7 (resp. 15) stages. An analysis of the local error coefficients supported by numerical 
tests shows that the methods proposed here are more efficient than those obtained in [CCDV09, HO09b] 
by applying the recursive triple jump and quadruple jump consti^uctions. An additional requirement when 
choosing a given method is that the arguments of the complex coefficients of the scheme have also to be 
taken into account. This constitutes a critical point for evolution equations where one of the operators (say, 
A ) has non-real eigenvalues in the right-hand side of the complex plane, as occurs, in particular, with the 
complex Ginzburg-Landau equation (1.6). In such a case, splitting methods of the form (1.2) where one of 
the two sets of coefficients 's or bi 's is entirely contained in the positive real axis, whereas the other set 
is included in the right-hand side of the complex planes are particularly well suited. Such splitting methods 
cannot be constructed as composition methods with the Strang splitting as basic method, so that a separate 
study is required to get the most efficient schemes within this class. 

Based on the theoretical framework worked out in [HO09a], the integrators proposed here can be applied 
to the numerical integration of linear evolution equations involving unbounded operators in an infinite dimen- 
sional space, like linear diffusion equations. As a matter of fact, although the theory developed in [HO09a] 
does not cover the generalization to semi-linear evolution equations, we have also included in our numerical 
tests a system of ODEs obtained from semi-linear evolution equations with a certain space discretization. 
All the numerical tests carried out with periodic boundary conditions show a considerable improvement in 
efficiency of our new methods with respect to existing splitting schemes. A remarkable feature of the new 
eighth-order composition method when applied to both the linear and semi-linear diffusion examples is that 
it is more efficient than all the other integrators of order p < 8 in the whole range of tolerances explored 
when periodic boundary conditions are considered. 

Concerning other (e.g. Dirichlet and Neumann) boundary conditions, the experiments carried out in 
[HO09b] for linear problems with methods obtained by applying the triple and quadruple jump technique 
show the existence of an order reduction phenomenon. Its origin is attributed by the authors of [HO09b] to 
the fact that Condition 2 in the introduction is not generally satisfied in this setting. In other words, terms 
of the form Ep^i e*("^+^) in (1.5) are not uniformly bounded on the interval [0, T] for some T > . As 
a consequence, the classical convergence order is no longer guaranteed in that case. This order reduction 
is also present, of course, when the methods introduced in this paper are applied to linear but also semi- 
linear problems. Nevertheless, as has been pointed out in [HO09b], splitting schemes of high order involving 
complex coefficients may still be advantageous in comparison with, say, Strang splitting when applied to 
linear parabolic problems with Dirichlet or Neumann boundary conditions, as they often lead to smaller 
errors and the order reduction is somehow confined to a neighborhood of the boundary. The situation, in 
our opinion, calls for a detailed study of the effect of boundary conditions other than the periodic case in 
the global efficiency of splitting methods with complex coefficients and the analysis of their applicability for 
more general parabolic problems than those considered in this paper. 
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